Predicting Theft Risk in San Francisco

knitr::opts_chunk$set(
    message = FALSE,
    warning = FALSE,
    results=FALSE,
    fig.align="center"
)

library(tidyverse)
library(tidycensus)
library(sf)
library(kableExtra)
library(patchwork)
library(ggplot2)
library(RSocrata)
library(viridis)
library(spatstat)
library(raster)
library(spdep)
library(FNN)
library(grid)
library(gridExtra)
library(knitr)

options(scipen=999)
options(tigris_class = "sf")
options(tigris_use_cache = TRUE)
source("functions.r")

p5 = plasma(5)
p2 = plasma(2)
crs=7131

1. Why analyze theft crimes in San Francisco?

1.1 Status

According to FBI UCR data, San Francisco is one of the most crime-ridden city in 2019, with a extremely high property crime rate of 5506.0 crimes out of 100,000 people, which is 2.61 times the U.S. average.

This high property crime rate was making residents in San Francisco feel unsafe¹ , and may prevent migrant population growth and potential population loss from dangerous areas.

read.csv("./data/crime_by_city_2019.csv")%>%
  mutate(City=gsub("(Metropolitan|\\d|Police Department)","",City),
         State=na_if(State,""),
         Population=as.numeric(gsub(",","",Population)),
         PropertyCrimeRate =as.numeric(gsub(",","",Property.crime))/Population*100000)%>%
  fill(State)%>%
  filter(Population>500000) %>% na.omit()%>%
  mutate(City=reorder(City,PropertyCrimeRate,mean),
         toHighlight=ifelse(City=="San Francisco","",0),
         toHighlight=na_if(toHighlight,0))%>%
  ggplot(aes(PropertyCrimeRate,City))+
  geom_col(aes(fill=toHighlight),show.legend=F)+
  geom_vline(aes(xintercept=2109.9),linetype=2)+
  annotate("text", x=2200, y=1,label="US AVG",size=3,angle=90,vjust=1,hjust=0)+
  scale_fill_manual(values=p5[3], guide="none",na.value="grey75")+
  labs(x="Property crimes per 100,000 population",
       y="Cities(population>500,000)",
       title="Property Crime Rate of US Major Cities",
       subtitle="FBI UCR 2019")+
  plotTheme()

Data source: FBI 2019, https://ucr.fbi.gov/crime-in-the-u.s/2019/crime-in-the-u.s.-2019/tables/table-8/

1.2 Background

There are several possible reasons for the high crime rate in San Francisco:

1. Historical reason

   Property crime rate in San Francisco hasn’t seen clear improvement in the past 10 years (2010~2020), which implies that the high crime rate is a long-term historical problem for SF.

2. Low arrest rate

   San Francisco ranks high in crime rate while it ranks low in arrests.² Especially the defunding of the police department causes the labour shortage, so the arrests rate cannot keep pace with the high crime rate.

3. Active crime

   The active crimes may also account for the high crime rate. There may be more thieves in number, more active thieves or both.

“When you see an activity like that go up that much, you are either talking more active burglars, or more burglars” – UC Berkeley School of Law professor and criminologist Franklin Zimring.

Thus, it’s certainly necessary to model geospatial risk and bring predictive policing into life. In this project, a geospatial risk model is built to “borrow” the experience from the past theft crime data and forecast latent risk of thefts. There is an important underlying hypothesis that crime risk is a function of exposure to a series of geospatial risk and protective factors, such as graffiti or abandoned buildings.

2. Theft Data in San Francisco

2.1 Theft point map

A map of your outcome of interest in point form, with some description of what, when, and why you think selection bias may be an issue.

A map of theft crime is plotted with the data from DataSF. Unlike fire or property deals, theft crimes may have a much higher selection bias, which will probably leads to the systematic racism of model results.

The following conditions may result in selection bias:

• If selective enforcement of thefts is the result of an officer’s preconceived or biased beliefs, then this “selection bias” will be contained into the crime outcome.

• If the victim chooses not to report the crime because of modest theft value, hard or cannot reach the police, then this “selection bias” will be in the outcome.

crime.2018 = st_read("./data/crime.2018.precleaned.geojson")%>%
  st_transform(crs = crs)
crime.2019 = st_read("./data/crime.2019.precleaned.geojson")%>%
  st_transform(crs = crs)
nbh = st_read("./data/nbh-sf.geojson")%>%
  dplyr::select(name)%>%
  st_transform(crs=crs)%>%
  rename(nbh.name=name)

crime.2018%>%
  ggplot()+
  geom_sf(size=0.3,fill="#dddddd",data=nbh)+
  geom_sf(size=0.3,alpha=0.8,shape=1,color=p5[3])+
  labs(title="Thefts Point Map",
       subtitle="San Francisco, 2018")+
  mapTheme()

As is shown in the figure above, thefts are concentrated in northwest San Francisco, where the downtown locates, and they are more sparse in the rest of the city. In addition, the nature of property thefts determine that the spatial distribution of these points is in grids, corresponding to the city road image.

Data source: DataSF, https://data.sfgov.org/Public-Safety/Police-Department-Incident-Reports-2018-to-Present/wg3w-h783

2.2 Fishnet Map

# city boundary
boundary = st_union(nbh)%>%st_buffer(1.5)

# make fishnet with cellsize 150m (around 500ft)
fishnet =
  st_make_grid(boundary,
               cellsize = 150, # meters 
               square = TRUE) %>%
  .[boundary] %>%
  st_sf() %>%
  mutate(uniqueID = rownames(.))

# add a value of 1 to each crime, sum them with aggregate
crime_net =
  dplyr::select(crime.2018) %>% 
  mutate(count = 1) %>% 
  aggregate(., fishnet, sum) %>%
  mutate(count = replace_na(count, 0),
         uniqueID = rownames(.),
         cvID = sample(round(nrow(fishnet) / 57), # to make 100 groups
                       size=nrow(fishnet), replace = TRUE))

ggplot() +
  geom_sf(data = crime_net, aes(fill = count), color = NA) +
  scale_fill_viridis(option="C") +
  labs(title = "Count of Thefts for the fishnet",
       subtitle="San Francisco, 2018") +
  mapTheme()

Instead of 500 meters, a 150-meter cell size is used to create the fishnet because of the smaller block size in San Francisco and denser crimes. Another reason of shrinking the cell size is that the city size is also relatively small.

Clustered spikes of high-theft-count can be observed, implying the underlying spatial process.

2.3 Risk Factors Map

Point Level: Aggregate

# 311 data
var.311 = c("Graffiti","Abandoned Vehicle","Damaged Property","Homeless Concerns","Streetlights","Street and Sidewalk Cleaning")

data.311 = st_read("./data/311.cleaned.geojson")%>%
  st_transform(crs)%>%
  filter(service_name%in%var.311)%>%
  mutate(Legend=service_name)%>%
  dplyr::select(Legend)

predictors.net.1 = data.311%>%
  st_join(., fishnet, join=st_within) %>%
  st_drop_geometry() %>%
  group_by(uniqueID,Legend) %>%
  summarize(count = n()) %>%
  full_join(fishnet, by = "uniqueID") %>%
  spread(Legend, count, fill=0) %>%
  st_sf() %>%
  dplyr::select(-`<NA>`) %>%
  na.omit() %>%
  ungroup()

“Graffiti”,“Abandoned Vehicle”,“Damaged Property”,“Homeless Concerns”,“Streetlights”,“Street and Sidewalk Cleaning” are chosen as exposure factors according to the broken window theory.

Data Source: 311

Point Level: NN distance

nn = function(df,k=3){
  return(nn_function(st_coordinates(st_centroid(predictors.net.1)), 
            st_coordinates(df), k=k))
}

# San Francisco Municipal Transportation stations
muniStops=read.csv("./data/muni.stops.csv")%>%
  st_as_sf(wkt="shape",crs=4326)%>%
  st_transform(crs)%>%
  mutate(Legend="muniStops")%>%
  dplyr::select(Legend)

# Bay Area Rapid Transit stations
bartStops=st_read("./data/bart.station.geojson")%>%
  st_transform(crs)%>%
  mutate(Legend="bartStops")%>%
  dplyr::select(Legend)

# get 5 nearest neighbor
predictors.net.2 = predictors.net.1 %>%
  mutate(muniStops.nn=nn(muniStops),
         bartStops.nn=nn(bartStops))

# distance to downtown
downtownPoint = filter(nbh, nbh.name == "Downtown / Union Square") %>%
  st_centroid()

predictors.net.2$downtownDistance =
  st_distance(st_centroid(predictors.net.2),
              downtownPoint) %>%
  as.numeric() 

# calculate nn for 311 vars
for (n in 1:(var.311%>%length())) {
  vname=var.311[n]%>%paste0(".nn")
  predictors.net.2 = predictors.net.2%>%
    mutate(!!vname := !! data.311%>%filter(Legend==var.311[n])%>%nn())
}

Public transportation like BART(Bay Area Rapid Transit) stops and MUNI(San Francisco Municipal Transportation) stations are taken into account. All predictors are transformed into the NN(Nearest Neighbor) distance in order to improve the overall quality of the model.

Data Source:

• MUNI Stops

• BART Stations

Polygons Level: Catagories

# police districts
police = st_read("./data/police.geojson")%>%
  st_transform(crs)%>%
  dplyr::select(company)%>%
  rename(police.ID=company)

# generaal zoning district (residential, commercial ...)
zoning.general = st_read("./data/zoning.geojson")%>%
  st_transform(crs)%>%
  group_by(gen)%>%
  summarise()%>%
  rename(zoning=gen)

# get acs data
acs.var = c(
  "B02001_001E", #pop
  "B02001_002E", #white alone
  "B06009_005E", #bachelor
  "B06012_002E", #poverty
  "B19013_001E", #MedHHInc
  "B25058_001E", #medRent
  "B23025_007", #unemployment num
  "B23025_001E" #over16
)
acs18 =
  get_acs(geography = "tract", variables = acs.var,
          year = 2018, state=06, county=075, geometry=T,output="wide",cache_table = T) %>%
  st_transform(crs)  %>% 
  rename(TotalPop = B02001_001E, Whites = B02001_002E,
         Bachelors = B06009_005E,TotalPoverty = B06012_002E,
         MedHHInc = B19013_001E, MedRent = B25058_001E,
         unemployNum = B23025_007E,o16=B23025_001E) %>%
  mutate(pctWhite = ifelse(TotalPop > 0, Whites / TotalPop,0),
         pctBachelors = ifelse(TotalPop > 0, (Bachelors / TotalPop),0),
         pctPoverty = ifelse(TotalPop > 0, TotalPoverty / TotalPop, 0),
         Majority_White = ifelse(pctWhite>median(pctWhite,na.rm=T),1,0),
         popDensity = TotalPop/st_area(geometry)%>%as.numeric(),
         unemploymentRate = unemployNum/o16,
         MedHHInc = replace_na(MedHHInc,0))%>%
  dplyr::select(popDensity,pctWhite,pctBachelors,pctPoverty,Majority_White,
                MedRent,MedHHInc,unemploymentRate)%>%filter(popDensity>0)

# adding Polygons
predictors.net.3 =
  st_centroid(predictors.net.2) %>%
  #st_join(police, by = "uniqueID") %>%
  st_join(zoning.general, by = "uniqueID") %>%
  st_join(nbh, by = "uniqueID") %>%
  st_join(acs18, by = "uniqueID") %>%
  st_drop_geometry() %>%
  left_join(dplyr::select(predictors.net.2, geometry, uniqueID)) %>%
  st_sf()

Census tract data of 2018 such as percentage bachelor, percentage white, population density, unemployment rate are included in the model to provide urban context.

Data Source:

• Neighborhoods

• Zoning Code

• ACS Tracts

• Police Districts

Visualize the risk factors by fishnet created above.

var.names = predictors.net.3%>%colnames()%>%setdiff(c("uniqueID","geometry"))
plots1=list()
for (var in var.names){
  is.discrete = predictors.net.3%>%pull(var)%>%is.character()
  p = ggplot(predictors.net.3)+
    geom_sf(aes_(fill=as.name(var)),size=0.001,show.legend = !is.discrete)+
    scale_fill_viridis(option = "C",discrete = is.discrete, name="")+
    labs(title=var)+
    mapTheme()
  plots1[[var]]=p
}
grid.arrange(grobs = plots1,top=textGrob("Risk Factors by Fishnet", gp=gpar(fontsize=28)))

2.4 Local Moran’s I

# make polygon to neighborhoods... 
crime.net.nb <- poly2nb(as_Spatial(crime_net), queen=TRUE)
# make and neighborhoods to list of weigths
crime.net.weights <- nb2listw(crime.net.nb, style="W", zero.policy=TRUE)

# local Morans'I
local_morans <- localmoran(crime_net$count, crime.net.weights, zero.policy=TRUE) %>% 
  as.data.frame()

# join local Moran's I results to fishnet
crime.net.localMorans <- 
  cbind(local_morans, as.data.frame(crime_net)) %>% 
  st_sf() %>%
  dplyr::select(Thefts_Count = count, 
                Local_Morans_I = Ii, 
                P_Value = `Pr(z != E(Ii))`) %>%
  mutate(Significant_Hotspots = ifelse(P_Value <= 0.001, 1, 0)) %>%
  gather(Variable, Value, -geometry)

vars.moran = unique(crime.net.localMorans$Variable)
plots2 = list()
for(i in vars.moran){
  plots2[[i]] = ggplot() +
      geom_sf(data = filter(crime.net.localMorans, Variable == i), 
              aes(fill = Value), colour=NA) +
      scale_fill_viridis(name="",option="C") +
      labs(title=i) +
      mapTheme(title_size = 14) + theme(legend.position="bottom")}
grid.arrange(grobs=plots2, ncol = 4, 
             top=textGrob("Local Morans I statistics", gp=gpar(fontsize=18)))

Figure above shows that relatively high values of Moran’s I, which represents strong and statistically significant evidence of spatial clustering of dependent variable. P-value and significant hot spot maps further shows the hot spot areas at the downtown area of SF.

2.5 Correlations scatter plots

final_net = left_join(crime_net, st_drop_geometry(predictors.net.3), by="uniqueID") 


final_net = final_net %>% 
  mutate(hotspot.isSig = ifelse(local_morans[,5] <= 0.001, 1, 0)) %>%
  mutate(hotspot.isSig.dist = 
           nn_function(st_coordinates(st_centroid(final_net)),
                       st_coordinates(st_centroid(filter(final_net, hotspot.isSig == 1))), k = 1))%>% 
  drop_na()


correlation.long= final_net%>%
  st_drop_geometry() %>%
    dplyr::select(-uniqueID, -cvID, -zoning,-nbh.name) %>%
    gather(Variable, Value, -count)

correlation.cor =
  correlation.long %>%
    group_by(Variable) %>%
    summarize(correlation = cor(Value, count, use = "complete.obs"))
    
ggplot(correlation.long, aes(Value, count)) +
  geom_point(size = 0.1,color=p5[4]) +
  geom_text(data = correlation.cor, aes(label = paste("r =", round(correlation, 2))),
            x=-Inf, y=Inf, vjust = 1.5, hjust = -.1) +
  geom_smooth(method = "lm", se = FALSE, colour = "black") +
  facet_wrap(~Variable, nrow = 4, scales = "free") +
  labs(title = "Theft count as a function of risk factors") +
  plotTheme(title_size = 28)

A small multiple scatterplot with correlations is shown above. Important factors are damaged property, damaged streetlights, clean requests and homeless concerns.

2.6 Histogram of theft count

ggplot(final_net) +
  geom_histogram(aes(count),bins=50,fill=p5[4],color="black") +
  labs(title = "Histogram of theft count",
       subtitle = "San Francisco, 2018") +
  plotTheme()

The histogram shows that the distribution of theft count per fishnet basically resembles Poisson distribution.

3. Regression

# define the variables we want
reg.ss.vars = "Abandoned Vehicle.nn,Damaged Property.nn,Graffiti.nn,Homeless Concerns.nn,Street and Sidewalk Cleaning.nn,Streetlights.nn,Street and Sidewalk Cleaning.nn,downtownDistance,popDensity,unemploymentRate,pctWhite,pctBachelors,MedRent,MedHHInc,hotspot.isSig,hotspot.isSig.dist"%>%str_split(",")%>%.[[1]]
reg.vars = "Abandoned Vehicle.nn,Damaged Property.nn,Graffiti.nn,Homeless Concerns.nn,Street and Sidewalk Cleaning.nn,Streetlights.nn,Street and Sidewalk Cleaning.nn,downtownDistance,popDensity,unemploymentRate,pctWhite,pctBachelors,MedRent,MedHHInc"%>%str_split(",")%>%.[[1]]
crossvalidation_id = "nbh.name"

# RUN REGRESSIONS
reg.cv = crossValidate(
    dataset = final_net,
    id = "cvID",                           
    dependentVariable = "count",
    indVariables = reg.vars) %>%
  dplyr::select(cvID , count, Prediction, geometry)

reg.ss.cv = crossValidate(
    dataset = final_net,
    id = "cvID",                           
    dependentVariable = "count",
    indVariables = reg.ss.vars) %>%
  dplyr::select(cvID , count, Prediction, geometry)

reg.spatialCV = crossValidate(
    dataset = final_net,
    id = crossvalidation_id,                           
    dependentVariable = "count",
    indVariables = reg.vars) %>%
  dplyr::select(cvID=as.name(crossvalidation_id) , count, Prediction, geometry)

reg.ss.spatialCV = crossValidate(
    dataset = final_net,
    id = crossvalidation_id,                           
    dependentVariable = "count",
    indVariables = reg.ss.vars) %>%
  dplyr::select(cvID=as.name(crossvalidation_id) , count, Prediction, geometry)

In the regression part, 4 different cross validation are used to compare the impact of factors on the final model. Specifically, two models are using spatial cross validation, i.e., to group the fishnet into neighborhoods, for each group, train the model on the other groups and test on selected model. The other two models are using groups randomly generated (about 100 groups). Besides, the variables in the models are different, spatial variables are controlled over the 4 cross validations.

4. Results

4.1 Multiple Map of Model Errors

# calculate errors by regression
reg.summary = rbind(
    mutate(reg.cv,           Error = Prediction - count,
                             Regression = "Random k-fold CV: Just Risk Factors"),
                             
    mutate(reg.ss.cv,        Error = Prediction - count,
                             Regression = "Random k-fold CV: Spatial Process"),
    
    mutate(reg.spatialCV,    Error = Prediction - count,
                             Regression = "Spatial LOGO-CV: Just Risk Factors"),
                             
    mutate(reg.ss.spatialCV, Error = Prediction - count,
                             Regression = "Spatial LOGO-CV: Spatial Process")) %>%
  st_sf()  

error_by_reg_and_fold =
  reg.summary%>%
  group_by(Regression,cvID) %>% 
  summarize(Mean_Error = mean(Prediction - count, na.rm = T),
            MAE = abs(Mean_Error)) %>%
  ungroup()

## plot histogram of OOF (out of fold) errors
error_by_reg_and_fold %>%
  ggplot(aes(MAE)) + 
    geom_histogram(bins = 30, colour="black", fill = p5[4]) +
    #scale_x_continuous(breaks = seq(0, 11, by = 1)) + 
    #xlim(c(0,11))+
    facet_wrap(~Regression) + 
    labs(title="Distribution of MAE", subtitle = "LOGO-CV",
         x="Mean Absolute Error", y="Count")+
    plotTheme()

As can be seen from the figure, there is little difference between the four CV methods. Among them, MAE of K-fold is mainly concentrated between 0 and 3 with fewer outliers, while MAE of LOGO has a higher frequency of 0, but there are more outliers, indicating that the prediction deviation of Spatial LOGO is very large for some areas, but smaller for others. The accuracy of LOGO is more extreme.

error_by_reg_and_fold %>%
  filter(str_detect(Regression, "LOGO")) %>%
  ggplot() +
    geom_sf(data=fishnet,fill=p5[1],size=0) +
    geom_sf(aes(fill = MAE)) +
    facet_wrap(~Regression) +
    scale_fill_viridis(option="C") +
    labs(title = "Theft errors by LOGO-CV Regression") +
    mapTheme() + theme(legend.position="bottom")

Just as the previous analysis, both models have large MAE in the downtown area, indicating that the models still have unexplained spatial processes. However, the spatial process cannot be explained by simply adding more spatial variables, and how to eliminate excessive MAE in downtown area is still a problem to be studied.

4.2 Table of MAE and standard deviation MAE by regression

k = st_drop_geometry(error_by_reg_and_fold) %>%
  group_by(Regression) %>% 
    summarize(Mean_MAE = round(mean(MAE), 2),
              SD_MAE = round(sd(MAE), 2)) %>%
  kable() %>%
    kable_styling("striped", full_width = F) %>%
    row_spec(2, color = "black", background = "#F89441aa") %>%
    row_spec(4, color = "black", background = "#F89441aa") 

k

Regression	Mean_MAE	SD_MAE
Random k-fold CV: Just Risk Factors	0.54	0.51
Random k-fold CV: Spatial Process	0.53	0.49
Spatial LOGO-CV: Just Risk Factors	1.21	2.84
Spatial LOGO-CV: Spatial Process	1.09	2.59

In general, random K-fold regression is better than LOGO regression, and spatial process plays a limited role in the model.

4.3 Raw errors by race context

t = reg.summary %>% 
  filter(str_detect(Regression, "LOGO")) %>%
    st_centroid() %>%
    st_join(acs18) %>%
    na.omit() %>%
      st_drop_geometry() %>%
      group_by(Regression, Majority_White) %>%
      summarize(mean.Error = mean(Error, na.rm = T)) %>%
      spread(Majority_White, mean.Error) %>%
      rename(Majority_non_White=`0`,Majority_White=`1`)%>%
      kable(caption = "Mean Error by tract racial context") %>%
        kable_styling("striped", full_width = F)  
t

Mean Error by tract racial context

Regression	Majority_non_White	Majority_White
Spatial LOGO-CV: Just Risk Factors	-0.0731254	0.0611867
Spatial LOGO-CV: Spatial Process	-0.0473177	0.0392849

In general, the prediction result of the model is larger than the actual value in the region of Majority Non White, and smaller than the actual value in the region of Majority White. Compared with common models, spatial model has smaller bias by race and better generalizability.

4.4 Kernel density and risk predictions for 2019

ppp = as.ppp(st_coordinates(crime.2018), W = st_bbox(final_net))
KD = density.ppp(ppp, 1000)

KDE_sf = as.data.frame(KD) %>%
  st_as_sf(coords = c("x", "y"), crs = st_crs(final_net)) %>%
  aggregate(., final_net, mean) %>%
  mutate(label = "Kernel Density",
         Risk_Category = ntile(value, 100),
         Risk_Category = case_when(
           Risk_Category >= 90 ~ "90% to 100%",
           Risk_Category >= 70 & Risk_Category <= 89 ~ "70% to 89%",
           Risk_Category >= 50 & Risk_Category <= 69 ~ "50% to 69%",
           Risk_Category >= 30 & Risk_Category <= 49 ~ "30% to 49%",
           Risk_Category >= 1 & Risk_Category <= 29 ~ "1% to 29%")) %>%
  cbind(
    aggregate(
      dplyr::select(crime.2018) %>% mutate(count = 1), ., sum) %>%
    mutate(count = replace_na(count, 0))) %>%
  dplyr::select(label, Risk_Category, count)

risk_sf =
  filter(reg.summary, Regression == "Spatial LOGO-CV: Spatial Process") %>%
  mutate(label = "Risk Predictions",
         Risk_Category = ntile(Prediction, 100),
         Risk_Category = case_when(
           Risk_Category >= 90 ~ "90% to 100%",
           Risk_Category >= 70 & Risk_Category <= 89 ~ "70% to 89%",
           Risk_Category >= 50 & Risk_Category <= 69 ~ "50% to 69%",
           Risk_Category >= 30 & Risk_Category <= 49 ~ "30% to 49%",
           Risk_Category >= 1 & Risk_Category <= 29 ~ "1% to 29%")) %>%
  cbind(
    aggregate(
      dplyr::select(crime.2018) %>% mutate(count = 1), ., sum) %>%
      mutate(count = replace_na(count, 0))) %>%
  dplyr::select(label,Risk_Category, count)

rbind(KDE_sf, risk_sf) %>%
  na.omit() %>%
  gather(Variable, Value, -label, -Risk_Category, -geometry) %>%
  ggplot() +
    geom_sf(aes(fill = Risk_Category), colour = NA) +
    geom_sf(data = sample_n(crime.2019, 3000), size = .3, colour = "#00000088") +
    facet_wrap(~label, ) +
    scale_fill_viridis(discrete = TRUE,option = "C") +
    labs(title="Comparison of Kernel Density and Risk Predictions",
         subtitle="2018 theft risk predictions; 2019 thefts") +
    mapTheme()

Although the prediction accuracy of the two models is similar in the downtown area, the Risk Predictions accuracy is significantly higher in other areas while the Kernel Density is very vague for other areas.

4.5 Bar plot of prediction results comparison

rbind(KDE_sf, risk_sf) %>%
  st_set_geometry(NULL) %>% na.omit() %>%
  gather(Variable, Value, -label, -Risk_Category) %>%
  group_by(label, Risk_Category) %>%
  summarize(count = sum(Value)) %>%
  ungroup() %>%
  group_by(label) %>%
  mutate(Rate_of_test_set_crimes = count / sum(count)) %>%
    ggplot(aes(Risk_Category,Rate_of_test_set_crimes)) +
      geom_bar(aes(fill=label), position="dodge", stat="identity") +
      scale_fill_manual(values = p5[c(4,3)]) +
      labs(title = "Risk prediction vs. Kernel density, 2018 Thefts") +
      plotTheme() + theme(axis.text.x = element_text(angle = 45, vjust = 0.5))

5. Conclusion

The model is satisfactory with the MAE of down to about 0.5, which means a high accuracy of result prediction. Moreover, the model has stronger generalization ability across different neighborhood contexts after adding spatial variables. As a result, the model may help control San Francisco’s high crime rate. However, the prediction accuracy of the model for some downtown areas is not high, and the reasons for this phenomenon and what kind of spatial process can be used to correct this error remain to be further explored.

However, although the model has relatively high promotion ability, selection bias is still inevitable. I think the model should not be put into use until we are sure whether the actual use of it will cause more serious racism.

1: Sanfrancisco.cbslocal.com. 2021. Fears Arise As Burglaries In San Francisco Soar. [online] Available at: https://sanfrancisco.cbslocal.com/2021/05/24/crime-and-punishment-fears-arise-as-burglaries-in-san-francisco-soar/ [Accessed 29 October 2021].